1. Предобработка данных

theme_set(theme_minimal())

birth_weight_orig <- readRDS("very_low_birthweight.RDS")

birth_weight <- birth_weight_orig %>% 
  select_if( ~ sum(is.na(.)) < 100) %>% 
  na.omit() %>% 
  mutate(across(c(race, inout, delivery, sex, dead, twn, vent, pneumo, pda, cld), ~ as.factor(.)))

2. Числовые переменные

birth_weight %>%
  select(where(is.numeric))%>%
  pivot_longer(everything())%>%
  ggplot(aes(x = value))+
  geom_density(fill = "lightblue")+
  facet_wrap(~name, scales = "free")

Удалим выбросы:

remove_outlier <- function(column) {
  Q <- quantile(column, probs=c(.25, .75))
  iqr <- IQR(column)
  column <- replace(column, column < (Q[1] - 1.5*iqr) | column > (Q[2]+1.5*iqr), NA)
  return(column)
}

bw_clean <- birth_weight %>%
  mutate(across(where(is.numeric), ~ remove_outlier(.x)))

bw_clean %>%
  select(where(is.numeric))%>%
  pivot_longer(everything())%>%
  ggplot(aes(x = value))+
  geom_density(fill = "lightblue", na.rm = T)+
  facet_wrap(~name, scales = "free")

Раскрасим переменные в зависимости от значения inout:

bw_clean %>%
  select(where(is.numeric) | inout)%>%
  pivot_longer(!inout)%>%
  ggplot(aes(x = value, fill = inout))+
  geom_density(alpha = 0.2, na.rm = T)+
  facet_wrap(~name, scales = "free")

3. lowph vs inout

bw_clean  %>% 
  ggplot(aes(y = lowph, fill = inout))+
  geom_boxplot()

Для сравнения средних значений переменной lowph в двух группах, объем которых не меньше 80, используем t-тест с поправкой Уэлча (функция t_test с настройками по умолчанию).

library(rstatix)
bw_clean %>% t_test(lowph ~ inout)%>%flextable::flextable()

.y.

group1

group2

n1

n2

statistic

df

p

lowph

born at Duke

transported

443

81

5.777689

108.8782

0.0000000729

Между группами наблюдается статистически значимое различие в уровнях pH, что может быть связано с разницей в выживаемости тех младенцев, которые родились в Дьюке, и тех, которые были транспортированы (у родившихся в Дьюке, возможно, выше как pH, так и выживаемость). Однако, чтобы проверить данное предположение, необходимо отдельно сравнить уровни выживаемости между группами, ведь различия могут быть вызваны и иными причинами.

4. Корреляция

library(corrplot)

bw_num <- bw_clean  %>% 
  select(where(is.numeric) & !birth & !year & !exit)%>% 
  na.omit()

bw_num%>% 
  cor()%>% 
  corrplot.mixed(order = 'AOE')

library(GGally)

bw_num%>% 
  ggpairs(progress = F)

5. Иерархическая кластеризация

library(pheatmap)
bw_num_scaled <- scale(bw_num)
bw_num_dist <- dist(bw_num_scaled)

res.hc <- hclust(bw_num_dist, method = "ward.D2")
plot(res.hc)
grp <- cutree(res.hc, k = 3)
rect.hclust(res.hc, k = 3, border = 2:5)

Я бы выделила в полученной иерархической структуре три больших кластера (обведены рамками)

6. Иерархическая кластеризация + heatmap

pheatmap(bw_num_scaled, 
         show_rownames = FALSE, 
         clustering_distance_rows = bw_num_dist,
         clustering_method = "ward.D2", 
         cutree_rows = 3,
         cutree_cols = length(colnames(bw_num_scaled)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")

Данный график демонстрирует нам структуру наших данных. Мы видим, что детей можно разделить на три примерно равных по объему кластера, однако они достаточно неоднородны, что может говорить о том, что на самом деле нужно выделять другое количество групп.

Если же говорить о переменных, то наиболее близкими друг к другу являются вес при рождении и гестационный возраст, а наиболее далеким от всех остальных переменных - количество дней, проведенных в больнице. Оставшиеся три переменные имеют почти что одинаковое положение в данной структуре. Это может объясняться тем, что вес и возраст положительно скоррелированы, причем коэффициент их корреляции наибольший из всех попарных корреляций. А количество дней, проведенных в больнице, отрицательно скоррелировано со всеми остальными переменными.

7. РСА

Проведем РСА анализ:

library(FactoMineR)
bw.pca <- prcomp(bw_num_scaled, scale = F)
summary(bw.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6
## Standard deviation     1.5954 1.0661 0.8823 0.8512 0.73334 0.52677
## Proportion of Variance 0.4242 0.1894 0.1298 0.1208 0.08963 0.04625
## Cumulative Proportion  0.4242 0.6136 0.7434 0.8641 0.95375 1.00000

Мы видим, что первые две компоненты суммарно объясняют больше половины дисперсии (61%), а чтобы объяснить не меньше 90%, нужно использовать 5 компонент.

fviz_eig(bw.pca, addlabels = T)

Наибольший вклад в первые две компоненты вносят вес и гестационный возраст (наиболее значимые для первой компоненты), а также количество тромбоцитов (наибольший влад во вторую компоненту)

fviz_pca_var(bw.pca, col.var = "contrib")

contrib_1 <- fviz_contrib(bw.pca, choice = "var", axes = 1) # 1
contrib_2 <- fviz_contrib(bw.pca, choice = "var", axes = 2) # 2

ggarrange(plotlist = list(contrib_1, 
                          contrib_2),
          ncol = 1)

8. biplot

Построим biplot для РСА:

library(ggbiplot)

bw <- bw_clean  %>% 
  na.omit() %>% 
  rowid_to_column("id")

bipl <- ggbiplot::ggbiplot(bw.pca, 
         scale=0, 
         groups = as.factor(bw$dead), 
         ellipse = T,
         alpha = 0.3) +
  theme_minimal()

bipl

9. Plotly

А также переведем его в “plotly”:

library(plotly)

bipl <-bipl+
  geom_point(data = bipl$data, aes(text = paste0("id: ", bw$id)), alpha=0)

ggplotly(bipl, tooltip = "text", width = 800)

10.Интерпретация

Про проведении РСА нам удалось перейти в новую систему координат, избавившись от значительной скоррелированности в исходных переменных. При этом первые два измерения в новой системе объясняют 61% дисперсии данных (достаточно высокий показатель), а наибольший вклад в эти две новые координаты вносят вес, возраст и количество тромбоцитов. Мы визуализировали в новой плоскости совместное распределение данных и переменных, которые разбиваются на две группы: умерших и выживших. Делать вывод об ассоциации смертности с исходными переменными в данном случае не вполне корректно, так как мы видим лишь проекции многомерных данных и исходных переменных на двумерную плоскость, а значит теряем значительную часть информации.

11. UMAP

Приведем наши данные к размерности 2 с помощью UMAP и визуализируем результат:

library(tidymodels)
library(embed)

umap_prep <- recipe(~., data = bw_num) %>% # "техническая" строка, нужная для работы фреймворка tidymodels
  step_normalize(all_predictors()) %>% # нормируем все колонки
  step_umap(all_predictors()) %>%  # проводим в UMAP
  prep() %>%  # "техническая" строка, нужная для работы фреймворка tidymodels. Мы выполняем все степы выше 
  juice() # Финальная строка - приводим результаты UMAP к стандартизированному датасету

plot_umap <- umap_prep %>%
  ggplot(aes(UMAP1, UMAP2)) +
  geom_point(aes(color = as.character(bw$dead)),
             alpha = 0.7, size = 2) +
  labs(title = "Исходный вариант", color = NULL) 

plot_umap

12. Изменение параметров UMAP

Можно заметить, что при изменение параметров программы паттерн изображения точек может кардинально измениться:

umap_prep_2 <- recipe(~., data = bw_num) %>% # "техническая" строка, нужная для работы фреймворка tidymodels
  step_normalize(all_predictors()) %>% # нормируем все колонки
  step_umap(all_predictors(),
            neighbors = 5,
            min_dist = 0.05) %>%  # проводим в UMAP
  prep() %>%  # "техническая" строка, нужная для работы фреймворка tidymodels. Мы выполняем все степы выше 
  juice() # Финальная строка - приводим результаты UMAP к стандартизированному датасету


umap_prep_2 %>%
  ggplot(aes(UMAP1, UMAP2))+
  geom_point(aes(color = as.character(bw$dead)),
             alpha = 0.7, size = 2) +
  labs(title = "Вариант с измененными параметрами", color = NULL)

13. Пермутация и UMAP

При пермутации 50% и 100% одной из переменных меняется изображение на плоскости UMAP:

bw_num_50 <-  bw_num %>%
  slice(1:244) %>%
  mutate(bwt = sample(bwt))%>%
  bind_rows(bw_num %>%
  slice(245:488))

umap_prep_3  <- recipe(~., data = bw_num_50) %>% # "техническая" строка, нужная для работы фреймворка tidymodels
  step_normalize(all_predictors()) %>% # нормируем все колонки
  step_umap(all_predictors()) %>%  # проводим в UMAP
  prep() %>%  # "техническая" строка, нужная для работы фреймворка tidymodels. Мы выполняем все степы выше 
  juice() # Финальная строка - приводим результаты UMAP к стандартизированному датасету

plot_umap_50 <- umap_prep_3 %>%
  ggplot(aes(UMAP1, UMAP2)) +
  geom_point(aes(color = as.character(bw$dead)),
             alpha = 0.7, size = 2) +
  labs(title = "50% пермутации", color = NULL) 


bw_num_100 <-  bw_num%>%
  mutate(bwt = sample(bwt))

umap_prep_3  <- recipe(~., data = bw_num_100) %>% # "техническая" строка, нужная для работы фреймворка tidymodels
  step_normalize(all_predictors()) %>% # нормируем все колонки
  step_umap(all_predictors()) %>%  # проводим в UMAP
  prep() %>%  # "техническая" строка, нужная для работы фреймворка tidymodels. Мы выполняем все степы выше 
  juice() # Финальная строка - приводим результаты UMAP к стандартизированному датасету

plot_umap_100 <-  umap_prep_3 %>%
  ggplot(aes(UMAP1, UMAP2)) +
  geom_point(aes(color = as.character(bw$dead)),
             alpha = 0.7, size = 2) +
  labs(title = "100% пермутации", color = NULL) 


ggarrange(plotlist = list(plot_umap, 
                          plot_umap_50, 
                          plot_umap_100),
          ncol = 3)

А при РСА анализе сокращается доля суммарной дисперсии, объяснененной несколькими первыми главными компонентами. Скорее всего это связано с тем, что вес вносил наибольший вклад в первую компоненту.

bw_50.pca <- prcomp(bw_num_50, scale = T)
summary(bw_50.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5     PC6
## Standard deviation     1.4825 1.0573 0.8875 0.8645 0.8401 0.66574
## Proportion of Variance 0.3663 0.1863 0.1313 0.1246 0.1176 0.07387
## Cumulative Proportion  0.3663 0.5526 0.6839 0.8085 0.9261 1.00000
bw_100.pca <- prcomp(bw_num_100, scale = T)
summary(bw_100.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5     PC6
## Standard deviation     1.4009 1.0747 0.9628 0.8767 0.8449 0.68766
## Proportion of Variance 0.3271 0.1925 0.1545 0.1281 0.1190 0.07881
## Cumulative Proportion  0.3271 0.5196 0.6741 0.8022 0.9212 1.00000

Визуализация точек и biplot также меняются. Мы можем заметить, как с увеличением доли пермутирования сокращается вклад переменной bwt в первые две компоненты:

bipl_50 <- ggbiplot::ggbiplot(bw_50.pca, 
         scale=0, 
         groups = as.factor(bw$dead), 
         ellipse = T,
         alpha = 0.3) +
  theme_minimal()+
  labs(title = "50% пермутации")

bipl_100 <- ggbiplot::ggbiplot(bw_100.pca, 
         scale=0, 
         groups = as.factor(bw$dead), 
         ellipse = T,
         alpha = 0.3) +
  theme_minimal()+
  labs(title = "100% пермутации")


bipl <- bipl +
  labs(title = "Исходный вариант") 

ggarrange(plotlist = list(bipl, 
                          bipl_50, 
                          bipl_100),
          ncol = 3,
          common.legend = TRUE)

14-15. Анализ чувствительности

Повторим пункты 4-6 для датасета, в котором сразу удалили все строки с пропущенными значениями. При таком подходе все колонки сохранены, но значений получается намного меньше.

bw_drop_na <- birth_weight_orig %>% 
  na.omit() %>% 
  mutate(across(c(race, inout, twn, magsulf, meth, toc, delivery, sex, dead, vent, pneumo, pda, cld, pvh, ivh, ipe), ~ as.factor(.)))%>%
  mutate(across(where(is.numeric), ~ remove_outlier(.x))) %>% 
  na.omit()

Можно заметить, что изменились коэффициенты корреляции: какие-то увеличились, какие-то уменьшились.

bw_drop_na_num <- bw_drop_na  %>% 
  select(where(is.numeric) & !birth & !year & !exit)%>% 
  na.omit()

bw_drop_na_num%>% 
  cor()%>% 
  corrplot.mixed(order = 'AOE')

Иерархическая кластеризация строк глобально изменилась не сильно:

bw_num_scaled <- scale(bw_drop_na_num)
bw_num_dist <- dist(bw_num_scaled)

res.hc <- hclust(bw_num_dist, method = "ward.D2")
plot(res.hc)
grp <- cutree(res.hc, k = 3)
rect.hclust(res.hc, k = 3, border = 2:5)

А вот структура соотношения переменных немного другая, например, lowph больше похож на bwt чем apg1 или pltct.

pheatmap(bw_num_scaled, 
         show_rownames = FALSE, 
         clustering_distance_rows = bw_num_dist,
         clustering_method = "ward.D2", 
         cutree_rows = 3,
         cutree_cols = length(colnames(bw_num_scaled)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")

В РСА сильно упала доля дисперсии, объясняемая первыми компонентами (например, первые две объясняют только 54%, хотя раньше был 61%):

bw.pca <- prcomp(bw_num_scaled, scale = F)
summary(bw.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.6022 1.0944 1.0267 0.8849 0.82135 0.67412 0.51869
## Proportion of Variance 0.3667 0.1711 0.1506 0.1119 0.09637 0.06492 0.03843
## Cumulative Proportion  0.3667 0.5378 0.6884 0.8003 0.89665 0.96157 1.00000
fviz_eig(bw.pca, addlabels = T)

А группы выживших и умерших сильнее расходятся в новых координата (возможно это связано с тем, что группы стали меньше):

bipl_drop_na <- ggbiplot::ggbiplot(bw.pca, 
         scale=0, 
         groups = as.factor(bw_drop_na$dead), 
         ellipse = T,
         alpha = 0.3) +
  theme_minimal()+
  labs(title = "Датасет без удаления столбцов")


ggarrange(plotlist = list(bipl, 
                          bipl_drop_na),
          ncol = 2,
          common.legend = TRUE)

umap_prep <- recipe(~., data = bw_drop_na_num) %>% # "техническая" строка, нужная для работы фреймворка tidymodels
  step_normalize(all_predictors()) %>% # нормируем все колонки
  step_umap(all_predictors()) %>%  # проводим в UMAP
  prep() %>%  # "техническая" строка, нужная для работы фреймворка tidymodels. Мы выполняем все степы выше 
  juice() # Финальная строка - приводим результаты UMAP к стандартизированному датасету

plot_umap_drop_na <- umap_prep %>%
  ggplot(aes(UMAP1, UMAP2)) +
  geom_point(aes(color = as.character(bw_drop_na$dead)),
             alpha = 0.7, size = 2) +
  labs(title = "Датасет без удаления столбцов", color = NULL) 

ggarrange(plotlist = list(plot_umap, 
                          plot_umap_drop_na),
          ncol = 2)

Теперь попробуем заменить пропущенные значения на средние:

bw_mean_na <- birth_weight_orig %>% 
  mutate(across(c(race, inout, twn, magsulf, meth, toc, delivery, sex, dead, vent, pneumo, pda, cld, pvh, ivh, ipe), ~ as.factor(.)))%>%
  select_if( ~ sum(is.na(.)) < 100) %>% 
  mutate(across(where(is.numeric), ~ as.double(.)))%>%
  mutate(across(where(is.numeric), ~ replace_na(., mean(., na.rm=TRUE))))%>%
  mutate(across(where(is.numeric), ~ remove_outlier(.x))) %>% 
  na.omit()
bw_mean_na_num <- bw_mean_na  %>% 
  select(where(is.numeric) & !birth & !year & !exit)%>% 
  na.omit()

bw_mean_na_num %>% 
  cor()%>% 
  corrplot.mixed(order = 'AOE')

Коэффициенты корреляции изменились не сильно

в иерархической структуре строк один из классов начал чуть сильнее выделяться на фоне двух остальных:

bw_num_scaled <- scale(bw_mean_na_num )
bw_num_dist <- dist(bw_num_scaled)

res.hc <- hclust(bw_num_dist, method = "ward.D2")
plot(res.hc)
grp <- cutree(res.hc, k = 3)
rect.hclust(res.hc, k = 3, border = 2:5)

Структура и иерархия переменных (столбцов) осталась прежней:

pheatmap(bw_num_scaled, 
         show_rownames = FALSE, 
         clustering_distance_rows = bw_num_dist,
         clustering_method = "ward.D2", 
         cutree_rows = 3,
         cutree_cols = length(colnames(bw_num_scaled)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")

В РСА опять же главные компоненты объясняют меньшую долю дисперсии, если сравнивать с исходным вариантом

bw.pca <- prcomp(bw_num_scaled, scale = F)
summary(bw.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6
## Standard deviation     1.5624 1.0673 0.8991 0.8592 0.76867 0.53144
## Proportion of Variance 0.4068 0.1898 0.1347 0.1231 0.09847 0.04707
## Cumulative Proportion  0.4068 0.5967 0.7314 0.8545 0.95293 1.00000
fviz_eig(bw.pca, addlabels = T)

biplot выглядит аналогично исходному варианту, только вторая компонента поменяла направление и размер единичного отрезка, поэтому график ортажен по горизонтали и точки лежат в промежутке [-3, 3], а не [-2, 2]

bipl_mean_na <- ggbiplot::ggbiplot(bw.pca, 
         scale=0, 
         groups = as.factor(bw_mean_na$dead), 
         ellipse = T,
         alpha = 0.3) +
  theme_minimal()+
  labs(title = "Датасет с импутацией данных")


ggarrange(plotlist = list(bipl, 
                          bipl_mean_na),
          ncol = 2,
          common.legend = TRUE)

Результаты UMAP не совпадают, хотя можно заметить похожие паттерны взаимного расположения точек:

umap_prep <- recipe(~., data = bw_mean_na_num) %>% # "техническая" строка, нужная для работы фреймворка tidymodels
  step_normalize(all_predictors()) %>% # нормируем все колонки
  step_umap(all_predictors()) %>%  # проводим в UMAP
  prep() %>%  # "техническая" строка, нужная для работы фреймворка tidymodels. Мы выполняем все степы выше 
  juice() # Финальная строка - приводим результаты UMAP к стандартизированному датасету

plot_umap_mean_na <- umap_prep %>%
  ggplot(aes(UMAP1, UMAP2)) +
  geom_point(aes(color = as.character(bw_mean_na$dead)),
             alpha = 0.7, size = 2)+
  labs(title = "Датасет с импутацией данных", color = NULL) 

ggarrange(plotlist = list(plot_umap, 
                          plot_umap_mean_na),
          ncol = 2)